      PROGRAM  SINK    (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT)
      COMMON/BLK1/DUM1(1),I(4),X(4),DX(4),RHO(4),CI(4),
     1 MCP( 4),NCP( 4),TTCP(10, 4),CPT(10, 4),NP( 4),NK( 4),TABP(5, 4),
     2TABT(10, 4),TABK(50, 4),CK(50, 4),CP(50, 4),CKO(50, 4),CIT( 4,10)
      COMMON/BLK2/DUM2(1),TO(100),T(100),TT(100),TTF(100),A(100),B(100),
     1 C(100),D(100),E(100),Z(100),NSTAS,NLYERS
      COMMON/BLK3/DUM3(1),QC,QCP,QR,QRP,QRR,HE,HEP,HWO,HW,LCHAN,QAERO,
     1 QAEROP,TF,TB,TP,TFP,TBP,TAU,DTAU,IPRINT,QRAT,PRAT,NBP
     3 ,NQRAT(10),MQRAT(10),QRATT(10,10),TAUQRAT(10,10),NPRAT(10),
     4 MPRAT(10),TAUPRAT(10,10),PRATT(10,10),QTAB(10),TAUT(10),PTAB(10)
     5,TTAU(10),WTCI(10),AREA(10),WTBP,WTSUM
      COMMON/BLK4/DUM4(1),TABX(10),TABZ(5),TABY(50),PL,CS,ISTAS,
     1 IT,SIGMA,ALPHA,ALPHO,EPS,EPSO,EPSB,EPSBO,MEPSG,NEPSG,TTEPSG(10),
     2 EPSGT(10),MH,NH,TTH(10),HT(10),HF,HFO,
     3 NTSP,MTS,NTS,TAUTS(20),TST(20)
      COMMON/BLK5/DUM5(1),MPT,NPT,TAUP(10),PT(10),MHF,NHF,TTHF(10),
     1 HFT(10),MTB,CST(10),
     1 NTB,TAUTB(10),TBT(10),MEPS,NEPS,TTEPS(10),EPST(10),MEPSB,NEPSB,
     2 TTEPSB(10),EPSBT(10),MTF,NTF,TAUTF(10),TFT(10),MAL,NAL,TTAL(10),
     3 ALT(10),MHW,NHW,TTHW(15),HWT(15),MQC,NQC,TAUQC(20),QCT(20),MQR,
     4 NQR,TAUQR(20),QRT(20),MHE,NHE,TAUHE(20),HET(20)
     5 ,NCHAN,TAUO,DTAU1,TCHAN,ERROR,PFREQ,ENDTAU,NQP,CF
      DIMENSION TITLE(12)
      NAMELIST/DATA1/TAUO,DTAU1,ENDTAU,PFREQ,ERROR,SIGMA,NLYERS,I,X,
     1 ISTAS,RHO,MCP,NCP,TTCP,CPT,MAL,NAL,TTAL,ALT,MEPS,NEPS,TTEPS,
     2 EPST,MH,NH,TTH,HT,MEPSG,NEPSG,TTEPSG,EPSGT,MHF,NHF,TTHF,HFT,
     3 MEPSB,NEPSB,TTEPSB,EPSBT,NP,TABP,TABT,NK,TABK,MPT,NPT,TAUP,PT,
     4 T,MHW,NHW,TTHW,HWT,MQC,NQC,TAUQC,QCT,MQR,NQR,TAUQR,QRT,MHE,NHE,
     5 TAUHE,HET,MTB,NTB,TAUTB,TBT,MTF,NTF,TAUTF,TFT,NTSP,MTS,NTS,
     6 TAUTS,TST,TCHAN,LCHAN,NCHAN,NQP,CF,NQRAT,MQRAT,TAUQRAT,QRATT,
     7 MPRAT,NPRAT,TAUPRAT,PRATT,CIT,CST,WTCI,AREA
C
C   PUT ZEROS IN COMMON STORAGE LOCATIONS
      DO 51 M=1,1017
   51 DUM1(M)=0.0
      DO 52 M=1,1003
   52 DUM2(M)=0.0
      DO 53 M=1,526
   53 DUM3(M)=0.0
      DO 54 M=1,166
   54 DUM4(M)=0.0
      DO 55 M=1,332
   55 DUM5(M)=0.0
  200 READ(5,101)TITLE
      IF(EOF,5)400,500
  400 STOP
  101 FORMAT(12A6)
  500 WRITE(6,102)TITLE
  102 FORMAT(12A6////)
      READ(5,DATA1)
      IPRINT=1
      WRITE(6,DATA1)
C
C   INITIALIZE INTERNAL QUANTITIES
      IF(CF.EQ.0.)CF=1.
      IT = 1
      TAU = TAUO
      TAUOO=TAUO
      NSTAS=1
      NBP=1
      NFREQ = 0
      WTSUM=0
      I(1)=I(1)-1
      DO 10 N=1,NLYERS
      DX(N)=X(N)/FLOAT(I(N))
   10 NSTAS=NSTAS+I(N)
      DO 11 N=1,NSTAS
      TO(N)=T(N)
   11 TT(N) = T(N)
C
C   STARTING POINT FOR NEW CALCULATION AFTER LAYER THICKNESS HAS BEEN
C     CHANGED
 210  DO 211 N=1,NSTAS
      T(N)=TO(N)
  211 TT(N)=T(N)
      DO 151 N=1,NLYERS
  151  CI(N)=CIT(N,NBP)
      CS=CST(NBP)
      TAU=TAUOO
      TAUO=TAUOO
      NFREQ=0
      IPRINT=1
      DTAU = DTAU1
      DTAUP=DTAU
      TAUPR=TAUO+PFREQ
      NN=NQRAT(NBP)
      MM=MQRAT(NBP)
      DO 13 N=1,NN
      QTAB(N)=QRATT(N,NBP)
   13 TAUT(N)=TAUQRAT(N,NBP)
      NM=NPRAT(NBP)
      MN=MPRAT(NBP)
      DO 14 N=1,NN
      PTAB(N)=PRATT(N,NBP)
   14 TTAU(N)=TAUPRAT(N,NBP)
      IF(LCHAN.EQ.0) GO TO 202
      DX(LCHAN)=X(LCHAN)/ I(LCHAN)
C
C   STARTING POINT FOR NEW TIME STEP CALCULATION
  202 TAU=TAU+DTAU
C
C   ADJUST TIME STEP TO INSURE PRINTING EXACTLY AT PFREQ
      IF(TAUPR-TAU)402,403,403
  402 DTAUP=DTAU
      DTAU=TAUPR-(TAU-DTAUP)
      TAU=TAU-DTAUP+DTAU
  403 DO 12 N=1,NSTAS
   12 T(N)=TT(N)
      CALL FTLUP (TAUO,QC,MQC,NQC,TAU QC,QCT)
      CALL FTLUP(TAUO,QRAT,MM,NN,TAUT,QTAB)
      QC=QC*QRAT
      CALL FTLUP (TAUO,HE,MHE,NHE,TAU HE,HET)
      CALL FTLUP (TAUO,QR,MQR,NQR,TAU QR,QRT)
      CALL FTLUP(TAUO,TB,MTB,NTB,TAUTB,TBT)
      CALL FTLUP (TAUO,TF ,MTF,NTF,TAUTF,TFT)
C
C   STARTING POINT FOR NEW CALCULATION WHEN TIME STEP HAS BEEN HALVED
  60  CALL FTLUP (TAU,QCP,MQC,NQC,TAU QC,QCT)
      CALL FTLUP(TAU,QRAT,MM,NN,TAUT,QTAB)
      QCP=QCP*QRAT
      CALL FTLUP (TAU,HEP,MHE,NHE,TAU HE,HET)
      CALL FTLUP (TAU,QRP,MQR,NQR,TAU QR,QRT)
      CALL FTLUP(TAU,PL,MPT,NPT,TAUP,PT)
      CALL FTLUP(TAU,PRAT,MN,NM,TTAU,PTAB)
      PL=PL*PRAT
      CALL FTLUP(TAU,TBP,MTB,NTB,TAUTB,TBT)
      CALL FTLUP (TAU,TFP,MTF,NTF,TAUTF,TFT)
C
C   STARTING POINT FOR NEW CALCULATION WHEN ERROR CRITERION HAS NOT BEEN
C     MET
  201 CALL FTLUP (TT(1),HW,MHW,NHW,TTHW,HWT)
      CALL FTLUP (TT(1),ALPHA,MAL,NAL,TTAL,ALT)
      CALL FTLUP (TT(1),EPS,MEPS,NEPS,TTEPS,EPST)
      CALL FTLUP(TT(NSTAS),EPSB,MEPSB,NEPSB,TTEPSB,EPSBT)
      CALL FTLUP (TT(NSTAS),HF,MHF,NHF,TTHF,HFT)
  203 IF(IT.NE.1)GO TO 207
      HWO=HW
      EPSO=EPS
      ALPHO=ALPHA
      EPSBO=EPSB
      HFO=HF
C
C   CALCULATE MATRIX COEFFICIENTS
  207 CALL COEFFS
C
C   SOLVE MATRIX
      CALL SOLMAT
C
C   TEST SOLUTION AGAINST ERROR CRITERION
   16 DO 45 N=1,NSTAS
      TEST=ABS((TT(N)-TTF(N))/TT(N))
      IF (TEST-ERROR) 45,45,17
  45  CONTINUE
      GO TO 23
  17  DO 19 N=1,NSTAS
      IF(IT.EQ.1)TT(N)=TTF(N)
   19 TT(N)=(TTF(N)+TT(N))/2.
      IT = IT+1
C
C   IF CALCULATION HAS ITERATED MORE THAN 9 TIMES, HALVE TIME STEP AND
C     TRY AGAIN
      IF(IT-9) 20,20,21
  21  TAU=TAU-DTAU
      DTAU=DTAU/2.
      TAU=TAU+DTAU
      IT=1
      DO 41 N=1,NSTAS
  41  TT(N)=T(N)
   20 GO TO 60
   23 DO 24 N=1,NSTAS
   24 TT(N)=TTF(N)
      TAUO=TAU
C
C   IF CALCULATION HAS ITERATED LESS THAN 3 TIMES, DOUBLE TIME STEP
      IF(IT-2)303,303,304
  303 DTAU=2.*DTAU
C
C   DO NOT LET TIME STEP GET TOO LARGE
      IF(DTAU.GT.10.) DTAU=10.
      IF(DTAU.GT.PFREQ) DTAU=PFREQ
  304 IT=1
  301 IF (NFREQ) 25,25,26
C
C      INITIAL PRINT AFTER FIRST TIME STEP
   25 CALL ZPRINT
      NFREQ=1
   26 IPRINT = 2
C
C   CHECK FOR TEMPERATURE LIMIT CRITERION
      IF (TCHAN-1.) 32,30,30
   30 TCK=.1*(ENDTAU-TAUOO)+TAUOO
C
C   IF THE TIME IS LESS THAN 10 PERCENT OF THE TOTAL TIME DO NOT CHECK
C     FOR PEAK IN BACK SURFACE TEMPERATURE
      IF(TAU.LT.TCK) GO TO 32
C
C   IF BACK SURFACE TEMPERATURE IS LESS THAN FRONT SURFACE TEMPERATURE
C     DO NOT CHECK FOR PEAK IN BACK SURFACE TEMPERATURE
      IF(T(1)-T(NSTAS))31,32,32
C
C   IF BACK SURFACE TEMPERATURE HAS PEAKED, PRINT
   31 IF(T(NSTAS)-TTF(NSTAS)) 32,33,33
  33  CALL ZPRINT
      IF(TCHAN.EQ.0.)GO TO 34
      TCT=T(NCHAN)/TCHAN
C
C   DETERMINE IF STATION TEMPERATURE IS WITHIN 1 PERCENT OF SPECIFIED
C     TEMPERATURE
      IF(.99.LT.TCT.AND.TCT.LT.1.01)34,35
C
C   IF TEMPERATURE LIMIT CRITERION IS MET, START CALCULATION FOR NEW
C     BODY POINT
  34  WTBP=0
      DO 150 N=1,NLYERS
      WT=X(N)*RHO(N)
  150 WTBP=WTBP+WT
      WTSUM=WTSUM+(WTBP+WTCI(NBP  ))*AREA(NBP  )
      CALL ZPRINT
      IF(NBP-NQP)300,200,200
  300 NBP=NBP+1
C
C   EXTRAPOLATED GUESS AT LAYER THICKNESS FOR NEW BODY POINT CALCULATION
      X(LCHAN)=X(LCHAN)*(QRATT(1,NBP  )/QRATT(1,NBP  -1))**.5
      GO TO 210
C
C   IF TEMPERATURE LIMIT CRITERION IS NOT MET, CHANGE LAYER THICKNESS
C     AND BEGIN NEW CALCULATION
  35  TRT=(T(NCHAN)-TCHAN)/(TCHAN-TO(NCHAN))
      IF(TRT.LT.-.5) TRT=-.5
      IF(TRT.GT.1) TRT=1
      X(LCHAN)=X(LCHAN)+X(LCHAN)*TRT*CF
      GO TO 210
  32  IF (ENDTAU-TAU) 33,33,27
   27 IF(ABS(TAUPR-TAU).LT.1.E-8)29,   202
   29 CALL ZPRINT
      TAUPR=TAUPR+PFREQ
      DTAU=DTAUP
      GO TO 202
   28 CALL ZPRINT
      GO TO 200
      END
      SUBROUTINE COEFFS
      COMMON/BLK1/DUM1(1),I(4),X(4),DX(4),RHO(4),CI(4),
     1 MCP( 4),NCP( 4),TTCP(10, 4),CPT(10, 4),NP( 4),NK( 4),TABP(5, 4),
     2TABT(10, 4),TABK(50, 4),CK(50, 4),CP(50, 4),CKO(50, 4),CIT( 4,10)
      COMMON/BLK2/DUM2(1),TO(100),T(100),TT(100),TTF(100),A(100),B(100),
     1 C(100),D(100),E(100),Z(100),NSTAS,NLYERS
      COMMON/BLK3/DUM3(1),QC,QCP,QR,QRP,QRR,HE,HEP,HWO,HW,LCHAN,QAERO,
     1 QAEROP,TF,TB,TP,TFP,TBP,TAU,DTAU,IPRINT,QRAT,PRAT,NBP
     3 ,NQRAT(10),MQRAT(10),QRATT(10,10),TAUQRAT(10,10),NPRAT(10),
     4 MPRAT(10),TAUPRAT(10,10),PRATT(10,10),QTAB(10),TAUT(10),PTAB(10)
     5,TTAU(10),WTCI(10),AREA(10),WTBP,WTSUM
      COMMON/BLK4/DUM4(1),TABX(10),TABZ(5),TABY(50),PL,CS,ISTAS,
     1 IT,SIGMA,ALPHA,ALPHO,EPS,EPSO,EPSB,EPSBO,MEPSG,NEPSG,TTEPSG(10),
     2 EPSGT(10),MH,NH,TTH(10),HT(10),HF,HFO,
     3 NTSP,MTS,NTS,TAUTS(20),TST(20)
      DIMENSION TCP(10),CCPT(10)
      IF(ISTAS.GT.0) I1STAS=ISTAS-1.
      L=1
      K=1
C
C   DETERMINE SPECIFIC HEAT AND THERMAL CONDUCTIVITY VALUES
      DO 7 N=1,NLYERS
      MMCP=MCP(N)
      NNCP=NCP(N)
      DO 200 M=1,NNCP
      TCP(M)=TTCP(M,N)
  200 CCPT(M)=CPT(M,N)
      NX=NK(N)/NP(N)
      NZ=NP(N)
      NY=NK(N)
      DO 310 KK=1,NY
  310 TABY(KK)=TABK(KK,N)
      DO 320 KK=1,NZ
  320 TABZ(KK)=TABP(KK,N)
      DO 330 KK=1,NX
  330 TABX(KK)=TABT(KK,N)
      K=K+I(N)
      DO 6 J=L,K
      TAV=(T(J)+TT(J))/2.
      CALL FTLUP (TAV,CCP,MMCP,NNCP,TCP,CCPT)
      CP(J,N)=CCP
      IF(J.EQ.NSTAS)TT(J+1)=TT(J)
      TAVG=(TT(J)+TT(J+1))/2.0
      IF(NZ.GT.1)GO TO 331
      CALL FTLUP(TAVG,CCK,1,NY,TABX,TABY)
      GO TO 332
  331 CALL DISCOT(TAVG,PL,TABX,TABY,TABZ,11,NY,NZ,CCK)
  332 CK(J,N)=CCK
      IF(IT.EQ.1)CKO(J,N)=CK(J,N)
    6 CONTINUE
    7 L=K
C
C   CALCULATE COEFFICIENTS AT FRONT SURFACE
   23 FAC1=(3.*CK(1,1))/(2.*DX(1)**2)
      FAC2=CK(2,1)/(6.*DX(1)**2)
      FAC4= RHO (1)*CP(1,1)+8./(3.*DX(1))*CS
      A(1)=0.
      B(1)=FAC4/DTAU+FAC1+4.*SIGMA*EPS*(TT(1)**3)/(3.*DX(1))
      C(1)=-FAC1-FAC2
      E(1)=FAC2
      DFAC2=(3.*CKO(1,1)*(T(2)-T(1)))/(2.*DX(1)**2)-(CKO(2,1)*(T(3)-T(2)
     1))/(6.*DX(1)**2)
      QRR=SIGMA*EPSO*T(1)**4
      QAEROP=QCP*(1.-HW/HEP)+ALPHA*QRP
      QAERO=QC*(1.-HWO/HE)+ALPHO*QR
      DFAC3=QAEROP+QAERO-QRR
      D(1)= FAC4*T(1)/DTAU+DFAC2+4./(3.*DX(1))*DFAC3
    9 LSTAS=1
      JJ=2
      DO 60 N=1,NLYERS
      LSTAS=LSTAS+I(N)
      DO 30 J=JJ,LSTAS
   15 IF(J.EQ.NSTAS)GO TO 50
      IF(J.EQ.I1STAS) GO TO 35
      IF(J.EQ.ISTAS) GO TO 36
      IF(J.EQ.LSTAS)GO TO 40
C
C   CALCULATE COEFFICIENTS FOR SOLID LAYER
      DEN=2.*DX(N)**2
      FAC4=RHO (N)*CP(J,N)
  42  FAC1=(CK(J,N)+CK(J-1,N))/DEN
      A(J)=-CK(J-1,N)/DEN
      B(J)=FAC4/DTAU+FAC1
      C(J)=-CK(J,N)/DEN
      E(J)=0.
      FAC2=(CKO(J,N)*(T(J+1)-T(J))-CKO(J-1,N)*(T(J)-T(J-1)))/DEN
      D(J)=FAC4*T(J)/DTAU+FAC2
C
C   CALCULATE COEFFICIENTS AT SOLID-SOLID INTERFACE
      GO TO 30
   40 DEN1=6.*DX(N)
      DEN2=2.*DX(N)
      DEN3=2.*DX(N+1)
      DEN4=6.*DX(N+1)
      FAC1=3.*CK(J-1,N)/DEN2
      FAC2=(RHO(N)*CP(J,N)*DX(N)+RHO(N+1)*CP(J,N+1)*DX(N+1)+8./3.*
     1 CI(N))
      FAC3=3.*CK(J,N+1)/DEN3
      Z(J)=CK(J-2,N)/DEN1
      A(J)=-FAC1-CK(J-2,N)/DEN1
      B(J)=FAC2/DTAU+FAC1+FAC3
      C(J)=-FAC3-CK(J+1,N+1)/DEN4
      E(J)=CK(J+1,N+1)/DEN4
      D(J)=FAC2*T(J)/DTAU-3.*CKO(J-1,N)/DEN2*(T(J)-T(J-1))+3.*CKO(J,N+1)
     1 /DEN3*(T(J+1)-T(J))+CKO(J-2,N)/DEN1*(T(J-1)-T(J-2))-CKO(J+1,N+1)
     2 /DEN4*(T(J+2)-T(J+1))
      GO TO 80
C
C   CALCULATE COEFFICIENTS AT FRONT SIDE OF AIR GAP
   35 DEN1=8.*DX(N)
      CALL FTLUP(TT(J),EGAP,MEPSG,NEPSG,TTEPSG,EPSGT)
      IF(IT.EQ.1) EGAPO=EGAP
      CALL FTLUP(TT(J),H,MH,NH,TTH,HT)
      IF(IT.EQ.1) HO=H
      FAC1=CK(J-2,N)/DEN1
      FAC2=9.*CK(J-1,N)/(8.*DX(N))
      FAC3=CK(J,N+1)/DX(N+1)
      FAC4=3.*DX(N)*RHO(N)*CP(J,N)/(4.*DTAU)+DX(N+1)*RHO(N+1)*CP(J,N+1)
     1 /DTAU+2.*CI(N)/DTAU
      Z(J)=FAC1
      A(J)=-FAC2-FAC1
      B(J)=FAC4+FAC2+FAC3+H+SIGMA*EGAP*TT(J)**3
      C(J)=-FAC3-H-SIGMA*EGAP*TT(J+1)**3
      D(J)=FAC4*T(J)     -9.*CKO(J-1,N)*(T(J)-T(J-1))/(8.*DX(N))
     1 +CKO(J-2,N)*(T(J-1)-T(J-2))/(8.*DX(N))+CKO(J,N+1)*(T(J+1)-T(J))/
     2 DX(N+1)-HO*(T(J)-T(J+1))-SIGMA*EGAPO*(T(J)**4-T(J+1)**4)
      GO TO 80
C
C   CALCULATE COEFFICIENTS AT BACK SIDE OF AIR GAP
   36 DEN=8.*DX(N+1)
      FAC1=CK(J-1,N)/DX(N)
      FAC2=9.*CK(J,N+1)/DEN
      FAC3=CK(J+1,N+1)/DEN
      FAC4=DX(N)*RHO(N)*CP(J,N)/DTAU+3.*DX(N+1)*RHO(N+1)*CP(J,N+1)/
     1 (4.*DTAU)+2.*CI(N)/DTAU
      A(J)=-FAC1-H-SIGMA*EGAP*TT(J-1)**3
      B(J)=FAC4+FAC1+FAC2+H+SIGMA*EGAP*TT(J)**3
      C(J)=-FAC2-FAC3
      E(J)=FAC3
      D(J)=FAC4*T(J)     -CKO(J-1,N)*(T(J)-T(J-1))/DX(N)+9.*CKO(J,N+1)*
     1 (T(J+1)-T(J))/DEN-CKO(J+1,N+1)*(T(J+2)-T(J+1))/DEN+HO*(T(J-1)
     2 -T(J))+SIGMA*EGAPO*(T(J-1)**4-T(J)**4)
C
C   CALCULATE COEFFICIENTS AT BACK SURFACE
      GO TO 30
   50 DEN1=2.*DX(N)*DX(N)
      DEN2=3.*DEN1
      DEN3=4.*DX(N)
      DEN4=3.*DX(N)
      FAC1=RHO(N)*CP(J,N)+8.*CI(N)/DEN4
      FAC2=CK(J-2,N)/DEN2
      FAC3=3.*CK(J-1,N)/DEN1
      Z(J)=FAC2
      A(J)=-FAC2-FAC3
      B(J)=FAC1/DTAU+FAC3+4./DEN4*(SIGMA*EPSB*TT(J)**3+HF)
      D(J)=FAC1*T(J)/DTAU+CKO(J-2,N)*(T(J-1)-T(J-2))/DEN2-3.*CKO(J-1,N)*
     1 (T(J)-T(J-1))/DEN1+4./DEN4*(SIGMA*(-EPSBO*(T(J)**4-TB**4)+EPSB*
     2 TBP**4)+HF*TFP-HFO*(T(J)-TF))
C
C   ELIMINATE COEFFICIENT BELOW TRI-DIAGONAL RESULTING FROM USE OF
C     3 STATION BACKWARD DIFFERENCE
   80 A(J)=A(J)-(B(J-1)*Z(J))/A(J-1)
      B(J)=B(J)-(C(J-1)*Z(J))/A(J-1)
      D(J)=D(J)-(D(J-1)*Z(J))/A(J-1)
   30 CONTINUE
   60 JJ=LSTAS+1
C
C   OPTION FOR SPECIFYING A STATION TEMPERATURE AS A FUNCTION OF TIME
      IF(NTSP) 21,21,22
   22 CALL FTLUP (TAU,TTF(NTSP),MTS,NTS,TAUTS,TST)
      Z(NTSP)=0.
      A(NTSP)=0.
      B(NTSP)=1.
      C(NTSP)=0.
      E(NTSP)=0.
      D(NTSP)=TTF(NTSP)
   21  RETURN
      END
      SUBROUTINE SOLMAT
      COMMON/BLK2/DUM2(1),TO(100),T(100),TT(100),TTF(100),A(100),B(100),
     1 C(100),D(100),E(100),Z(100),NSTAS,NLYERS
      DIMENSION SB(100),W(100),G(100),DC(100),SC(100)
    8 SB(1) = C(1)/B(1)
      W(1) = B(1)
      G(1) = D(1)/W(1)
      SC(1)=E(1)/W(1)
      DO 4 N=2,NSTAS
      W(N) = B(N)-A(N)*SB(N-1)
      SC(N)=E(N)/W(N)
      SB(N)=(C(N)-A(N)*SC(N-1))/W(N)
    4 G(N) = (D(N)-A(N)*G(N-1))/W(N)
   14 N=NSTAS
      TTF(N)=G(N)
      N=NSTAS-1
      TTF(N)=G(N)-SB(N)*TTF(N+1)
    6 N=N-1
      TTF(N)=G(N)-SB(N)*TTF(N+1)-SC(N)*TTF(N+2)
      IF(N-1)7,7,6
    7 RETURN
      END
      SUBROUTINE ZPRINT
      COMMON/BLK1/DUM1(1),I(4),X(4),DX(4),RHO(4),CI(4),
     1 MCP( 4),NCP( 4),TTCP(10, 4),CPT(10, 4),NP( 4),NK( 4),TABP(5, 4),
     2TABT(10, 4),TABK(50, 4),CK(50, 4),CP(50, 4),CKO(50, 4),CIT( 4,10)
      COMMON/BLK2/DUM2(1),TO(100),T(100),TT(100),TTF(100),A(100),B(100),
     1 C(100),D(100),E(100),Z(100),NSTAS,NLYERS
      COMMON/BLK3/DUM3(1),QC,QCP,QR,QRP,QRR,HE,HEP,HWO,HW,LCHAN,QAERO,
     1 QAEROP,TF,TB,TP,TFP,TBP,TAU,DTAU,IPRINT,QRAT,PRAT,NBP
     3 ,NQRAT(10),MQRAT(10),QRATT(10,10),TAUQRAT(10,10),NPRAT(10),
     4 MPRAT(10),TAUPRAT(10,10),PRATT(10,10),QTAB(10),TAUT(10),PTAB(10)
     5,TTAU(10),WTCI(10),AREA(10),WTBP,WTSUM
      GO TO (20,30),IPRINT
   20 WRITE(6,126)
  126 FORMAT(1H1)
   30 WRITE(6,127)TAU,DTAU
  127 FORMAT(1H0,2X5HTAU =F12.6,10X6HDTAU =F20.8/2X17HTEMPERATURE TABLE)
      LN=1
      LL=1
      DO 10 N=1,NLYERS
      LL=I(N)+LL
      WRITE(6,129)(TT(J),J=LN,LL)
      LN=LL
   10 CONTINUE
  129 FORMAT(6E20.8)
      WRITE(6,128)QC,QCP,QR,QRP,QRR,HE,HEP,HWO,HW,X(LCHAN),QAERO,QAEROP,
     1 TF,TB,QRAT,PRAT,WTCI(NBP),AREA(NBP),WTBP,WTSUM
  128 FORMAT(1H0,1X6HQC   =E15.8,5X6HQCP  =E15.8,5X6HQR   =E15.8,5X6HQRP
     1  =E15.8,5X6HQRR  =E15.8/2X6HHE   =E15.8,5X6HHEP  =E15.8,5X6HHWO
     2=E15.8,5X6HHW   =E15.8,5X6HX    =E15.8/2X6HQAERO=E15.8,4X7HQAEROP=
     3E15.8,5X6HTF   =E15.8,5X6HTB   =E15.8,5X6HQRAT =E15.8/2X6HPRAT =
     4E15.8,5X6HWTCI =E15.8,5X6HAREA =E15.8,5X6HWTBP =E15.8,5X6HWTSUM=
     5E15.8)
 9999 RETURN
      END
      SUBROUTINE FTLUP (X,Y,M,N,VARI,VARD)                              E1.1  10
*********DOCUMENT DATE  7/7/69     SUBROUTINE REVISED  7/7/69  *********E1.1  20
*        MODIFICATION OF LIBRARY INTERPOLATION SUBROUTINE  FTLUP        E1.1  30
      DIMENSION VARI(1),VARD(1),V(3),YY(2)                              E1.1  40
      DIMENSION II(43)                                                  E1.1  50
*                                                                       E1.1  60
*      INITIALIZE ALL INTERVAL POINTERS TO -1.0   FOR MONOTONICITY CHECKE1.1  70
      DATA (II(J),J=1,43)/43*-1/                                        E1.1  80
      MA=IABS(M)                                                        E1.1  90
*                                                                       E1.1 100
*            ASSIGN INTERVAL POINTER FOR GIVEN VARI TABLE               E1.1 110
*      THE SAME POINTER WILL BE USED ON A GIVEN VARI TABLE EVERY TIME   E1.1 120
      LI=MOD(LOCF(VARI(1)),43)+1                                        E1.1 130
      I=II(LI)                                                          E1.1 140
      IF (I.GE.0) GO TO 10                                              E1.1 150
      IF (N.LT.2) GO TO 10                                              E1.1 155
*                                                                       E1.1 160
* MONOTONICITY CHECK                                                    E1.1 170
      IF (VARI(2)-VARI(1)) 1,1,3                                        E1.1 180
*  ERROR IN MONOTONICITY                                                E1.1 190
    2 K=LOCF (VARI(1))                                                  E1.1 200
      PRINT 102,J,K,(VARI(J),J=1,N),(VARD(J),J=1,N)                     E1.1 210
  102 FORMAT (1H1,* TABLE BELOW OUT OF ORDER FOR FTLUP  AT POSITION  *  E1.1 220
     1,I5,/* X TABLE IS STORED IN LOCATION *,O6,//(8G15.8))             E1.1 230
      STOP                                                              E1.1 240
*  MONOTONIC DECREASING                                                 E1.1 250
    1 DO 5 J=2,N                                                        E1.1 260
      IF (VARI(J)-VARI(J-1))5,2,2                                       E1.1 270
    5 CONTINUE                                                          E1.1 280
      GO TO 10                                                          E1.1 290
*  MONOTONIC INCREASING                                                 E1.1 300
    3 DO 6 J=2,N                                                        E1.1 310
      IF (VARI(J)-VARI(J-1))2,2,6                                       E1.1 320
    6 CONTINUE                                                          E1.1 330
*                                                                       E1.1 340
* INTERPOLATION                                                         E1.1 350
   10 IF (I.LE.0) I=1                                                   E1.1 360
      IF (I.GE.N) I=N-1                                                 E1.1 365
      IF (N.LE.1) GO TO 8                                               E1.1 370
      IF (MA.NE.0) GO TO 99                                             E1.1 380
*  ZERO ORDER                                                           E1.1 390
    8 Y=VARD(1)                                                         E1.1 400
      GO TO 800
*  LOCATE I INTERVAL (X(I).LE.X.LT.X(I+1))                              E1.1 430
   99 IF ((VARI(I)-X)*(VARI(I+1)-X)) 61,61,40                           E1.1 440
*  IN GIVES DIRECTION FOR SEARCH OF INTERVALS                           E1.1 450
   40 IN=SIGN(1.0,(VARI(I+1)-VARI(I))*(X-VARI(I)))                      E1.1 460
*  IF X OUTSIDE ENDPOINTS, EXTRAPOLATE FROM END INTERVAL                E1.1 470
   41 IF ((I+IN).LE.0) GO TO 61
      IF ((I+IN).GE.N) GO TO 61
      I=I+IN                                                            E1.1 490
      IF ((VARI(I)-X)*(VARI(I+1)-X)) 61,61,41                           E1.1 500
   61 IF (MA.EQ.2) GO TO 200                                            E1.1 510
*                                                                       E1.1 520
* FIRST ORDER                                                           E1.1 530
      Y=(VARD(I)*(VARI(I+1)-X)-VARD(I+1)*(VARI(I)-X))/(VARI(I+1)-VARI(I)E1.1 540
     1   )                                                              E1.1 550
      GO TO 800
*                                                                       E1.1 580
* SECOND ORDER                                                          E1.1 590
  200 IF (N.EQ.2) GO TO 2                                               E1.1 600
      IF (I.EQ.(N-1)) GO TO 209                                         E1.1 605
      IF (I.EQ.1) GO TO 201                                             E1.1 610
*  PICK THIRD POINT                                                     E1.1 620
      SK= VARI(I+1)-VARI(I)                                             E1.1 630
      IF ((SK*(X-VARI(I-1))).LT.(SK*(VARI(I+2)-X))) GO TO 209           E1.1 640
  201 L=I                                                               E1.1 650
      GO TO 702                                                         E1.1 660
  209 L=I-1                                                             E1.1 670
  702 V(1)=VARI(L)-X                                                    E1.1 680
      V(2)=VARI(L+1)-X                                                  E1.1 690
      V(3)=VARI(L+2)-X                                                  E1.1 700
      YY(1)=(VARD(L)*V(2)-VARD(L+1)*V(1))/(VARI(L+1)-VARI(L))           E1.1 710
      YY(2)=(VARD(L+1)*V(3)-VARD(L+2)*V(2))/(VARI(L+2)-VARI(L+1))       E1.1 720
      Y=(YY(1)*V(3)-YY(2)*V(1))/(VARI(L+2)-VARI(L))                     E1.1 730
  800 II(LI)=I
      RETURN                                                            E1.1 750
      END                                                               E1.1 760
      SUBROUTINE DISCOT (XA,ZA,TABX,TABY,TABZ,NC,NY,NZ,ANS)             E1.2   1
********* DOCUMENT DATE 08-01-68   SUBROUTINE REVISED 08-01-68 *********E1.2   2
C     THE DIMENSIONS IN THIS SUBROUTINE ARE ONLY DUMMY DIMENSIONS.      E1.2   3
      DIMENSION TABX(2),TABY(2),TABZ(2),NPX(8),NPY(8),YY(8)             E1.2   4
C     DIMENSION TABX(2),TABY(2),TABZ(2),NPX(8),NPY(8),YY(8)             E1.2   5
      CALL UNS (NC,IA,IDX,IDZ,IMS)                                      E1.2   6
      IF (NZ-1)   5,5,10                                                E1.2   7
    5 CALL DISSER (XA,TABX(1),1,NY,IDX,NN)                              E1.2   8
      NNN=IDX+1                                                         E1.2   9
      CALL LAGRAN (XA,TABX(NN),TABY(NN),NNN,ANS)                        E1.2  10
      GOTO 70                                                           E1.2  11
   10 ZARG=ZA                                                           E1.2  12
      IP1X=IDX+1                                                        E1.2  13
      IP1Z=IDZ+1                                                        E1.2  14
      IF (IA)   15,25,15                                                E1.2  15
   15 IF (ZARG-TABZ(NZ))   25,25,20                                     E1.2  16
   20 ZARG=TABZ(NZ)                                                     E1.2  17
   25 CALL DISSER (ZARG,TABZ(1),1,NZ,IDZ,NPZ)                           E1.2  18
      NX=NY/NZ                                                          E1.2  19
      NPZL=NPZ+IDZ                                                      E1.2  20
      I=1                                                               E1.2  21
      IF (IMS)   30,30,40                                               E1.2  22
   30 CALL DISSER (XA,TABX(1),1,NX,IDX,NPX(1))                          E1.2  23
      DO 35 JJ=NPZ,NPZL                                                 E1.2  24
      NPY(I)=(JJ-1)*NX+NPX(1)                                           E1.2  25
      NPX(I)=NPX(1)                                                     E1.2  26
   35 I=I+1                                                             E1.2  27
      GOTO 50                                                           E1.2  28
   40 DO 45 JJ=NPZ,NPZL                                                 E1.2  29
      IS=(JJ-1)*NX+1                                                    E1.2  30
      CALL DISSER (XA,TABX(1),IS,NX,IDX,NPX(I))                         E1.2  31
      NPY(I)=NPX(I)                                                     E1.2  32
   45 I=I+1                                                             E1.2  33
   50 DO 55 LL=1,IP1Z                                                   E1.2  34
      NLOC=NPX(LL)                                                      E1.2  35
      NLOCY=NPY(LL)                                                     E1.2  36
   55 CALL LAGRAN(XA,TABX(NLOC),TABY(NLOCY),IP1X,YY(LL))                E1.2  37
      CALL LAGRAN (ZARG,TABZ(NPZ),YY(1),IP1Z,ANS)                       E1.2  38
   70 RETURN                                                            E1.2  39
      END                                                               E1.2  40
      SUBROUTINE UNS (IC,IA,IDX,IDZ,IMS)                                E1.2  41
      IF (IC)   5,5,10                                                  E1.2  42
    5 IMS=1                                                             E1.2  43
      NC=-IC                                                            E1.2  44
      GOTO 15                                                           E1.2  45
   10 IMS=0                                                             E1.2  46
      NC=IC                                                             E1.2  47
   15 IF (NC-100)   20,25,25                                            E1.2  48
   20 IA=0                                                              E1.2  49
      GOTO 30                                                           E1.2  50
   25 IA=1                                                              E1.2  51
      NC=NC-100                                                         E1.2  52
   30 IDX=NC/10                                                         E1.2  53
      IDZ=NC-IDX*10                                                     E1.2  54
      RETURN                                                            E1.2  55
      END                                                               E1.2  56
      SUBROUTINE DISSER (XA,TAB,I,NX,ID,NPX)                            E1.2  57
      DIMENSION TAB(2)                                                  E1.2  58
C     DIMENSION TAB(2)                                                  E1.2  59
      NPT=ID+1                                                          E1.2  60
      NPB=NPT/2                                                         E1.2  61
      NPU=NPT-NPB                                                       E1.2  62
      IF (NX-NPT)   10,5,10                                             E1.2  63
    5 NPX=I                                                             E1.2  64
      RETURN                                                            E1.2  65
   10 NLOW=I+NPB                                                        E1.2  66
      NUPP=I+NX-(NPU+1)                                                 E1.2  67
      DO 15 II=NLOW,NUPP                                                E1.2  68
      NLOC=II                                                           E1.2  69
      IF (TAB(II)-XA)   15,20,20                                        E1.2  70
   15 CONTINUE                                                          E1.2  71
      NPX=NUPP-NPB+1                                                    E1.2  72
      RETURN                                                            E1.2  73
   20 NL=NLOC-NPB                                                       E1.2  74
      NU=NL+ID                                                          E1.2  75
      DO 25 JJ=NL,NU                                                    E1.2  76
      NDIS=JJ                                                           E1.2  77
      IF (TAB(JJ)-TAB(JJ+1))   25,30,25                                 E1.2  78
   25 CONTINUE                                                          E1.2  79
      NPX=NL                                                            E1.2  80
      RETURN                                                            E1.2  81
   30 IF (TAB(NDIS)-XA)   40,35,35                                      E1.2  82
   35 NPX=NDIS-ID                                                       E1.2  83
      RETURN                                                            E1.2  84
   40 NPX=NDIS+1                                                        E1.2  85
      RETURN                                                            E1.2  86
      END                                                               E1.2  87
      SUBROUTINE LAGRAN (XA,X,Y,N,ANS)                                  E1.2  88
      DIMENSION X(2),Y(2)                                               E1.2  89
C     DIMENSION X(2),Y(2)                                               E1.2  90
      SUM=0.0                                                           E1.2  91
      DO 3 I=1,N                                                        E1.2  92
      PROD=Y(I)                                                         E1.2  93
      DO 2 J=1,N                                                        E1.2  94
      A=X(I)-X(J)                                                       E1.2  95
      IF (A)   1,2,1                                                    E1.2  96
    1 B=(XA-X(J))/A                                                     E1.2  97
      PROD=PROD*B                                                       E1.2  98
    2 CONTINUE                                                          E1.2  99
    3 SUM=SUM+PROD                                                      E1.2 100
      ANS=SUM                                                           E1.2 101
      RETURN                                                            E1.2 102
      END                                                               E1.2 103
1SAMPLE CASE
 $DATA1
 TAUO=0,  DTAU1=.03125,  ENDTAU=1000,  PFREQ=200,ERROR=.001,
 EPSGT=.8,  NCHAN=23,       NQRAT=1,1,
 1,4,5,6,  MQRAT=0,0,0,1,1,1,  QRATT(1,1)=1.,  QRATT(1,2)=.8,  QRATT(1,3)=.6,
 NQP=6,
 QRATT(1,4)=.4,.4,.5,.4,  QRATT(1,5)=.2,.2,.25,.3,.2,  QRATT(1,6)=.1,.1,.12,.15,
 .15,.1,  TAUQRAT(1,4)=0,194.,200.,1000.,  TAUQRAT(1,5)=0,183.,190.,200.,1000.,
 TAUQRAT(1,6)=0,172.,180.,190.,250.,1000.,
 PRATT(1,1)=1.,   PRATT(1,2)=.64,  PRATT(1,3)=.36,  PRATT(1,4)=.16,
 PRATT(1,5)=.04,  PRATT(1,6)=.01,
 EPST=.8,     MCP=1,1,1,  NCP=2,1,1,  NP=3,1,1,  NK=6,1,1,  TABP=.01,.1,1.,
 MPT=1,  NPT=4,  PT=.06,.025,.1,1.,
 NLYERS=3,  I=20,1,2,  ISTAS=21,
 MHW=1,   NHW=10,   MQC=1,   NQC=15,
          TAUQC=0,26,45,61,66,74,88,92,103,114,142,170,197,229,1000.,
 MHE=1,   NHE=15,
          TAUHE=0,26,45,61,66,74,88,92,103,114,142,170,197,229,1000.,
 AREA=.1,.15,.2,.2,.2,.15,
 LCHAN=1,          TAUP=0.,70.,220.,1000.,       TAUQRAT(1,4)=0.,114.,120.,1000,
 TAUQRAT(1,5)=0.,103.,110.,120.,1000.,   TAUQRAT(1,6)=0.,92.,100.,110.,180.,
 1000.,
         SIGMA=5.66961E-11, X=.03048,.3048E-02,.5090E-03,  RHO=9.6,.8,2723.1,
 TTCP=275,1670,    CPT(1,1)=.6276,.12552E+1,  CPT(1,2)=1.046,  CPT(1,3)=.9205,
 TABT=275,1670.,  TABK(1,1)=.0312E-03,.2496E-03,.04368E-03,.3744E-03,.04992E-03,
 .4368E-03,  TABK(1,2)=.0312E-03,  TABK(1,3)=.23088,  T=40*294.,
 TTHW=278,555.,833.,1111.,1389.,1667.,1944.,2222.,2500.,2778.,  HWT=267.,548.,
 844.,1160.,1489.,1833.,2188.,2564.,3011.,3269.,  QCT=11.35,34.05,79.45,181.6,
 261.05,340.5,567.5,567.5,454.,340.5,227.,113.5,22.7,11.35,4.54,
 HET=397.,629.,954.,1366.,1580.,1868.,2460.,2460.,2153.,1868.,1357.,937.,617.,
 395.,290.,   TCHAN=422.,  CIT(1,1)=.41,  CIT(1,2)=.41,  CIT(1,3)=.41,
 CIT(1,4)=.41,  CIT(1,5)=.41,  CIT(1,6)=.41,
                                       CST=2.,2.,1.63,1.63,1.22,1.22,
 RHO=96.,
 WTCI=4.88,4.88,3.9,3.9,2.93,2.93,
 $
